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ABSTRACT 

The shearing sheet is a model dynamical system that is used to study the small-scale dy- 
namics of astrophysical disks. Numerical simulations of particle trajectories in the shearing 
sheet usually employ the leapfrog integrator, but this integrator performs poorly because of 
velocity-dependent (Coriolis) forces. We describe two new integrators for this purpose; both 
are symplectic, time-reversible and second-order accurate, and can easily be generalized to 
higher orders. Moreover, both integrators are exact when there are no small-scale forces such 
as mutual gravitational forces between disk particles. In numerical experiments these integra- 
tors have errors that are often several orders of magnitude smaller than competing methods. 
The first of our new integrators (SEX) is well-suited for disks in which the typical inter-particle 
separation is large compared to the particles' Hill radii (e.g., planetary rings), and the second 
(SEKI) is designed for disks in which the particles are on bound orbits or the separation is 
smaller than the Hill radius (e.g., irregular satellites of the giant planets). 

Key words: methods: N-body simulations; methods: numerical; celestial mechanics; planets: 
rings; planets and satellites: formation; 



2 HILL'S APPROXIMATION 

For simplicity we consider mainly the three-body problem, al- 
though the results we describe are easy to generahze to other sys- 
tems such as planetary rings (see Sect.[5]l. Thus we follow the mo- 
tion of two nearby small bodies with masses mi and m2 in the grav- 
itational field of a large body of mass M » oti, OT2 (a more careful 
version of this derivation is given by |Henon & Pe"titl|l986) . The 
two small bodies follow approximately the same orbit around the 
large body, and we assume that this mean orbit is circular with 
semi-major axis a. The angular speed of the mean orbit is then 
fl = [G(M+m)/a']''- where m = mi+m2- Formally, to derive Hill's 
equations of motion, one lets m/M shrink to zero while assuming 
that the mass ratio milmj is fixed and the separation between the 
two small bodies is of order (m/My^^. 

Let us adopt a local right-handed Cartesian coordinate system 
with its origin at mi, rotating uniformly with angular velocity SI. 
The unit vector e, points away from the larger mass M, the unit 
vector e,, points in the direction of motion along the mean orbit, 
and the unit vector e. points in the direction of the angular velocity 
vector of the mean orbit. Hill's equations of motion for the relative 
position r = r2 - Ti are then 

r = -2Q. e.xr + 3Q.^{r ■ e.,-) 6;, - fl-(r • e.) e- -i- f , (1) 

where the force exerted by nii on m2 is f = -VO with <t>(r) = 
Gm/\r\. The Hamiltonian coiTesponding to these equations of mo- 



1 INTRODUCTION 

Hill's approximation, or the shearing sheet approximation, is an 
essential tool for the study of the small-scale dynamics of astro- 
physical disks (for a general review of Hill's approximation see, 
e.g., |Binn ey & Tremaine 20081. The method was originally de- 
vised by [Hill (1878 1 to study the motion of the Moon, and has 
been applied to galaxy di sks (e.g.,[Go ldreich & Lynden-Bell 1965 
[Julian & Toomre| |1966', 'Goldreich & Tremaine 1978^, accretion 



disks (e.g., Hawley & Balbus 1992; Stone & Gardine r|2010|>, plan- 



etary rings (e.g.. Wisdom & Tremaine 1988;Salo 1992; Richardson 
1 1 994[ [Cnda et al. 20 1 , R ein & Papaloizo u 2010) and planetesimal 
disks (e.g.,|Tanga et ai.|2 004; Johansen et al.|2009[ |Bai & Stone] 
|2010[|Rein et al.|2010^ . Hill's approximation is widely used in nu- 
merical simulations when it is impossible to model an entire disk 
with adequate numerical resolution. 

This paper discusses numerical methods for following orbits 
in Hill's approximation. We first describe the relevant equations of 
motion in Sect. |2] We then review existing integrators in Sect. |3] 
and describe two new algorithms. In Sect. |4] we compare the con- 
vergence and performance of all these integrators. We summarize 
and discuss several generalizations in Sect.|5] 
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tion IS 

H(r, p) = + (p X r) • + [r' - 3 (r • e,)'] + 0(r) 

= /fo(r,p) + 0(r), (2) 

where p = r- flrxe-is the canonical momentum conjugate to r. 

Note that when f = the equations of motion are invariant 
under the shear transformation 



(3) 



for arbitrary constants c,. and c,,. This invariance allows the use of 
periodic shearing boundary conditions for the study of the local 
dynamics of disks and is one reason why Hill's approximation is so 
useful. 

Alternatively, we can work in a non-rotating but accelerated 
reference frame with origin at mi. If the relative position in this 
frame is R = {X, Y, Z) and the (x, y, z) and (X, ¥, Z) reference frames 
coincide at ? = 0, then 



R = f22([r-e,(0]e,(0-R) + f, 
with 6^(0 = cos(f2f)ex + sin(nf)ey- The Hamiltonian is 



(4) 



H{R, P, = ^P^ + T^^^ {R" - 3[R • e.vCO]^} + ®(R) (5) 

where <I)(R) = 0(r) = -Gm/\r\ since |R| = |r|. 

The physical interpretation of Eqs. ^ and l|2j, and the opti- 
mum choice of integrator, depend on the relative size of the terms 
on the right side. This can be parametrized using the Hill radius, 
rm\ = (\Gin/Q}y^^ , which sets the relevant scale in the problem. 

• If the mass of the small bodies is negligible (m — » 0, rnin — > 0) 
their relative motion is simply epicyclic motion, which can be 
solved exactly (see e.g. |Henon & Petit|1986[|Binney & Tremaine| 
|2008[ or below). If m is sufficiently small, the motion can be re- 
garded as a perturbed epicycle. Quantitatively, this requires that 
one or more of the first three terms on the right side of Eq. ^ is 
much larger than |f| = Gm/r^, or 



''Hill 



» min [l,(£lrHm/v)"^], 



(6) 



where v is the velocity. One example is a planetary ring with 
semi-major axis a much less than the Roche limit ajt = 
2.46[3M/(4;rp,,)]''''; here M is the planet mass, mi and mi are 
the masses of the ring particles, and pp is the ring-particle den- 
sity. For two particles the separation between them cannot be less 
than twice their radius, which implies that r/rnin > » 1. 

If the force due to the mutual gravity of the small masses is large, 
we can interpret the solution to Eq. |TJ as a Keplerian orbit per- 
turbed by Coriolis and tidal forces (the terms proportional to Si 
and Or respectively). For a circular orbit of semi-major axis r 
the ratio of these perturbing forces to the Kepler force Gm/r 
is (r/fHiii)''^ and (r/rnin)^ respectively. One example is the ir- 
regular satellites of the giant planets of the solar system, which 
typically have r/rm\ - 0.1-0.5. 



[Binney & Tremaine||2008| . However, as we shall show below, 
leapfrog does not work well in Hill's approximation because there 
are velocity-dependent forces. Leapfrog is best suited for integrat- 
ing equations of motion of the form 



r = -V<I)(r) 



(V) 



where 0(r) is a (in general time-dependent) potential. Note that 
Eq. ^ can not be written in this form. The Hamiltonian corre- 
sponding to Eq. ^ is 

Hir, p)^^y+ 0(r) = //K.n(p) + 1>(r). (8) 

We assume that we have the position and velocit3{]of a particle at 
time /", r" = r(f") and v" = r(/") = p(/")- The goal is to approximate 
the new position and velocity at time = /" + At. Leapfrog is 
usually written as a chain of three operators, labeled Kick, Drift, 
Kick, applied successivel)]^ 

y«+i/2 = V' + ^Af f" Kick 



I," +1/2 



+ Afv"+'/- 
+ |Af f"+ 



Drift 

Kick, 



where f" = -V<l)(r"). The kick operator corresponds to follow- 
ing the trajectory exactly under the influence of the Hamiltonian 
0(r), the potential energy, and the drift operator corresponds to the 
Hamiltonian /fKin(p), the kinetic energy. Thus if we denote the op- 
erator for the exact evolution of a trajectory for a time-step At un- 
der an arbitrary Hamiltonian H by H(At), a single leapfrog step 
is the operator Hi(^]„(l At)^(At)Hi^i„(}j At). Each of these operators 
is symplectic since they are governed by a Hamiltonian, which im- 
plies that the leapfrog operator is also symplectic jSaha & Tremaine] 
|1992[ l. Moreover leapfrog is time-reversible and second-order, that 
is, the error after a single time-step is O(Ar'). 

When leapfrog is used in Hill's approximation, to integrate 
the equations of motion governed by Eq. ^ rather than Eq. (j7]l, 
the additional terms on the right side are added in the kick step, 
since they change the velocity. Thus we have (partly in component 
notation, for clarity): 



,1+1/2 



11+ 1/2 



v'; -I- i Af (3£i2 x" + 2n v;; + /'•) 
v;. + ^At[-2Qvi + f^:) 



+ iA?(-fi"z" -I- 



f") 



r" + At v"+'/^ 



vf "2 + |Af (3f2-x"+' + 2Qv"/"^ + f;*') 



vr"^ + iA?(-2nvr'" + /;"+') 



Kick 



Drift 



Kick 



^1/2+ iAf(-nV+'-Hy:"+'). 



However, because of the velocity-dependent Coriolis force on the 
right side, in this form leapfrog loses most of its desirable prop- 
erties: (i) it is neither symplectic nor time-reversible; (ii) it is only 
first-order accurate rather than second-order, that is, the error after a 
single time-step is 0(At^y, (iii) it leads to secular drifts in quantities 
that should be conserved. 

In principle these difficulties could be avoided by using 



3 INTEGRATORS 
3.1 Leapfrog integrator 

Many N-body simulations, in Hill's approximation and other con- 
texts, use the standard leapfrog integrator (e.g., |Springel||2005| 



For tlie Hamiltonian in Eq. |8J the canonical momentum p is the same as 
the time derivative of the position (the velocity). This is not the case for the 
Hamiitonian in Hiif's approximation, Eq. J2j. 

2 One can also apply the operators in a mfferent order: Drift, Kick, Drift. 
This does not change any conciusion in this paper. 
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leapfrog to integrate the equations of motion in the form ([4]l, which 
have no velocity-dependent forces so leapfrog is symplectic, time- 
reversible, and second-order. We do not pursue this option because 
(i) these equations are not invariant under the shear transformation, 
Eq. l[3]l, which limits their usefulness for the study of disks and 
rings; (ii) the method described in Sect. |3.5| below is better. 

There are schemes for constructing symplectic integrators for 
arbitrary Hamiltonians, but these are generally implicit and often 
algebraically complicated, at least for high-order schemes. The first 
explicit symplectic integrator for trajectories in Hill's approxima- 
tion is due to |Heggie| ( [2001[ l. However, Heggie's integrator is not 
time-reversible and is only first-order; we will not discuss it further 
because our numerical experiments show that it is not competitive 
with the integrators described below. 



3.2 Modified leapfrog integrator 

One can improve the standard leapfrog algorithm by using a pre- 
dictor step to approximate the velocity-dependent forces at the end 
of the time-step. This leads to the following scheme 



11+1/2 



,n+l/2 

y 

"+1/2 



v;: + i A/ (-2^ v'; + y;") Kick 



-H+l 

/l+l 

H+l 

y 

,>"+' 



v;: + Af(-2nv!;. + ,/;") 

r" + ^t v"+'/2 



"+1/2 



H+l/2 



I Af(-2nv!;. + ,/;"+') 
+ ^A?(-f22j;«+i + /:'-+'). 



Drift 



Kick 



Here v"^' is the predicted value of the velocity at time t + Af. This 
integrator is second-order, one order higher than standard leapfrog, 
but is neither time-reversible nor symplectic. This method is im- 
plemented in the particle codes Gasoline and PkdGRAV l jWadsley| 
let al.|2004| l. 

3.3 Symmetrized leapfrog 



[Mikkola & Merritt| ( |2006[ l describe a simple algorithm that converts 
any one-step integrator to a time-reversible integrator. We have ap- 
plied the Mikkola-Merritt algorithm to the standard (first-order) 
leapfrog integrator of Sect. |3.1[ thereby upgrading it to a time- 
reversible and second-order (but not symplectic) integrator. Time- 
reversibility endows integrators with most of the same desirable 
properties as symplecticity. 

The tests described below show that symmetrized leapfrog is 
not competitive with some of the other integrators discussed here. 
However, this elegant integrator could be profitably applied to sys- 
tems described by more complicated time-reversible Hamiltonians 
and also systems that are not governed by a Hamiltonian at all. 



al. integrator is symplectic, time-reversible, and accurate to second 
order. We refer the reader to the original paper for a derivation. 
Here, we merely list the final algorithm, which can also be written 
as three operators. Kick, Drift, KiclJ^ 

^/j+i/4 ^ v-;.- iAf(nV-/;) 



P'l. = V;. + IQx" + ^ At f^! 

n+l/2 _ ,,>i+l/4 



/1+1/2 



II+3/4 



y 

, ,«+ 1 



- fijc" - n (x" + Af vf 

Vl + ^^^ti-Sl^z" +f") 



r" + Afv"+'/^ 



n+3/4 



At (n2x"+' 



Kick 



Drift 



Kick 



3.5 Symplectic epicycle integrator (SEI) 

Mixed variable symplectic (MVS) schemes such as the Wisdom- 
Holman integrator have become the method of choice for long- 
term integrations of planetary orbits ( [Wisdom & Holman|1991[[Ki-| 
|noshita et aL|199l| . Like leapfrog, MVS schemes split the Hamil- 
tonian into two parts, H(r, p) = //^(r, p) + ^B(r, p), each of which 
is anal5'tically integrable. First, the trajectory is advanced under the 
influence of Hb for half a time-step, then Ha for a full time-step, 
then Hb again. In contrast to leapfrog, where H^ and Hb are the ki- 
netic and potential energy (cf. Eq. [Sj, in MVS schemes H^ and 
Hb are chosen so that \Hb\ \Ha\- Thus in the planetary case 
Ha is chosen to be the Kepler Hamiltonian while Hb represents 
the small gravitational forces from other planets. MVS integrators 
are symplectic (since the trajectory is advanced by a sequence of 
Hamiltonian maps) and time-reversible, and the eiTor per time-step 
is 0(At^)0{HB/HA). 

Interestingly, it is even easier to derive an MVS integrator for 
Hill's equations of motion than for Keplerian motion. As we show 
below, it is possible to solve for the epicyclic motion, that is the 
motion that is governed by the Hamiltonian Ho(r, p) in Eq. ([T|, in 
closed form with almost no computational effort (see also |Henon| 
|& Petit|[l986) . We use this result to construct two new MVS in- 
tegrators, one each for the two cases given at the end of Sect. [2] 
systems in which the forces f due to the gravity of the small bodies 
or other sources are relatively small (this section), and systems in 
which the gravity between a single pair of bodies dominates their 
motion (Sect.[3]6}. 

We first solve the equations of motion governed by Ho(r,p). 
To do that we shift the particle to a coordinate system where the 
particle's center of epicyclic motion is at the origin. We then do 
a rotation to account for the evolution in the epicycle. Finally we 
shift back and account for the shear. 

Using the same notation as abov^ the center of epicyclic mo- 



3.4 Quinn et al. integrator 

Recently, iQuinn et al.|pOTol l described an integrator that exhibits 3 -pj^g Q^ft step is the same as in standard leapfrog, but the Kick step is 

the desirable features of leapfrog despite the presence of velocity- modified. 

dependent forces in Hill's approximation. In particular, the Quinn et * Note that v is the velocity, not the canonical momentum. 
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tion of a particle is at 

We can tlien define 

=n(x"-x2) 
y: = {Q.(f-yl). 



(9) 



(10) 



Tlie evolution of these quantities during one time-step Af can be 
written as a rotation around the origin with an angle QAf, 



= x'; cos(n AO + >■;' sin(n AO 
= -xll sin(fi AO + cos(n AO- 



(11) 



Now we have only to undo the previous shift to the center of 
epicyclic motion and account for the shear to get the position and 
velocity at the new time f + Af: 



(12) 



v;^' =-2xr'-|xSn. 

The integration of the vertical motion can also be described by a 
rotation, so that the new vertical position and velocity at time f + Af 
are given by 



= z" cos(n AO + \-"-Or^ sin(f2 Af) 
= -z"f2 sin(f2 AO + v'l cos(n AO- 



(13) 



In some uses of Hill's approximation such as galactic disks, the 
epicycle and vertical frequencies may differ from the azimuthal 
frequency £1, but this generalization is easy to incorporate. The op- 
erator corresponding to the steps l[9)-(|13^ may be written in our 
notation as i?o(Af). 

Note that no function evaluation had to be performed dur- 
ing the entire step (i.e., there is no call to sqrt()). The sines and 
cosines appearing in the above equations are constant and the same 
for all particles. They can be pre-calculated at the beginning of the 
time-step or even at the beginning of the simulation if the time-step 
is fixed. All other operations are additions and multiplications. No 
significant additional storage is needed when there are many parti- 
cles. Also note that the integrator can be completely described by 
translations and rotations, making it an attractive choice for pro- 
grams running on graphic processors (GPUs). 

In long, high-accuracy integrations, round-off^ errors in the ro- 
tations in steps l[TTJ and l[T3j can cause a linear drift in energy and 
other integrals of motion, at arate |Aii/ii| ~ £-(f/Af), where £ is the 
machine precision, typically 2^^^ for double-precision arithmetic. 
An elegant solution to this problem is described in Appendix [A] 

An MVS integrator that includes additional forces due to a 
potential 0(r) may then be written as 



ffsEi(Ao = m\ht) ©(AO m\ht\ 



(14) 



where as usual <I>(Af) represents the kick step 

v"+' = v"-Af V(l)(r"+"2). (15) 
This integrator is symplectic, time-reversible, and second-order. 



and in contrast to the other integrators we have discussed so far, 
becomes exact as V<5> — > 0. More precisely, if the gravitational po- 
tential C) is 0(e) then the error of the SEI integrator after a single 
time-step is 0{e Af'), while the error of the Quinn et al. integrator 
is 0(Af3). 

As pointed out by [Quinn et al.|p010^ , numerical codes that 
implement collision detection usually assume that particles move 
along straight lines. In that case collision detection can be done ex- 
actly (although it is often done approximately). In contrast to the 
leapfrog and Quinn et al. integrators, the trajectory of a particle 
in SEI is not a straight line between kick steps. This might make 
collision detection harder. However, the curved trajectories are a 
real feature of the physics in Hill's approximation. Therefore, it 
does not make sense to choose an integrator that solves the equa- 
tions of motion incorrectly just to search for collisions along those 
incorrect trajectories exactly - it is better to detect collisions ap- 
proximately along exact trajectories than the reverse. Developing 
an efficient collision algorithm for curved trajectories of this kind 
is a research problem that needs further work. An obvious first step 
is to re-use the already implemented collision detection algorithms 
by approximating the trajectory as the line that joins the initial and 
final positions defined by the curved trajectory. This should work 
reasonably well so long as flAf <g: 1 . 



3.6 Symplectic epicycle-Kepler integrator (SEKI) 

The integrator described in the previous subsection is designed for 
the case where the forces due to the potential O are small compared 
to the forces that govern the epicyclic motion. We now describe 
an integrator for a situation in which the force due to the Kepler 
potential cD(r) = -Gm/ |r| is comparable to or stronger than the 
forces governing the epicyclic motion. 

We first note that one can integrate motion in the Kepler 
Hamiltonian //Kcp(r, p) = //Kin(p) + ©(r) = ^p^ - Gm/|r| ex- 
actly up to machine precision. Efficient methods for doing so are 
described by Wisdom & Holmm] jl991^ . Also note that one can 
rewrite Eq. ^ as 



H{r, p) = //o(r, p) + Ht,,,(r, p) - H^,,(p). 



(16) 



This motivates the following scheme, which we call symplectic 
epicycle-Kepler integrator (SEKI): 



^SEKl(AO 



(17) 



/fo(^AO i?Ki„(-|AO i?Kep(AO /?Ki„(-|AO Hoi^.At). 

Note that the drift operator, Hni„, has a negative time-step. This 
scheme is symplectic, second-order, and time-reversible. These 
statements are also true for other symmetric permutations of these 
operators. This particular permutation has been chosen because the 
computationally most expensive operator //kcp is called only once 
per time-step. Of course, if output is not needed at every time-step 
the half-steps at the end of step n and the start of step f! + 1 can be 
combined; then the alternative scheme 

/fKep(^AO i^K,„(-^AO Ho(At) /fKi„(-lAO iJKep(fAO (18) 

is hardly more expensive. 

Note that these operators are defined in (r, p) phase space. For 
the operators HKin and i?Kcp the canonical momentum p is equal to 
the velocity v = r = dH/dp but for Ho we have v = p + f2r x e-. 
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Figure 1. Pailicle motion during one epicycle using different integrators 
and a time-step = 0.1 ■ 2nQr^ . The initial position, which is also the 
exact final position, is marked by a circle. 



4 TESTS 

We ran many test integrations to study the convergence, accuracy, 
and computational cost of the integrators presented in the previous 
section. We present only four representative examples. All of these 
tests use the algorithm in Appendix|A]to minimize round-off errors. 
Without loss of generality, we set f2 = 1 from now on. 



4.2 Perturbed epicyclic motion 

The motion of a test particle in the presence of a mass m + 
can be described by a perturbed epicycle when the mass is suf- 
ficiently small, or, in other words, when the test particle is suf- 
ficiently far away from the mass or moving sufficiently fast, as 
quantified by Eq. |6](. We place the particle initially on a circular 
orbit around mass M (uniform motion along e,, in Hill's approxi- 
mation), with positions and velocities r = (5.55,2613.91,0)'^ and 
r = (0,-8.32,0)^. The trajectory passes the perturbing mass at 
an impact parameter corresponding to 8 rnin. The integration time 
is 100 epicycle periods. As an astrophysical example of such tra- 
jectories, we refer the reader to a study of the stochastic motion of 
moonlets embedded in Saturn's rings by |Rein & Papaloizou| ( |2010l l. 

In Figure [3] we plot the same diagnostics as in Figure [2] and 
additionally the phase error. We also looked at other measures of 
the accuracy of the integrators but do not show the results. SEI (for 
unbound orbits) and SEKI (for bound orbits, see below) perform at 
least as well as the other integrators tested by all measures that we 
have examined, and much better by most of them. SEI and SEKI 
also turn out to be exceptionally good at calculating the phase of 
the epicyclic motion. 

One can see from Figure [3^ that SEI and SEKI exhibit en- 
ergy errors that are up to three orders of magnitude smaller than 
any other integrator for typical time-steps used in most simulations 
(A/ ~ 10"''~5 • 10"^). Eventually the errors are dominated by round- 
off errors (A? < 5 ■ lO"'*) for all integrators. One can further see that 
at a fixed time-step SEI produces a phase error that is up to seven ( !) 
orders of magnitude smaller than the error produced by the Quinn 
et al. integrator. From FigurejSj) it is clear that SEI is also by far the 
fastest integrator for a given precision. 



4.1 Epicyclic motion 

We first examine the case in which the perturbing force f in Eq. |TJ 
vanishes, which corresponds to m ^ or rnin — > 0. In this case 
Hill's equations of motion can be solved analytically and lead to 
epicyclic motion. We initialize the particle position to r = (1, 0, 0)^ 
and the velocity to r = (0, -2,0)^^, which corresponds to a tra- 
jectory that is a clockwise closed ellipse centered on the origin. 
We integrate the trajectory forward in time for one epicycle period 
(/ = 2;rn-'). 

To illustrate the behavior of each integrator, we use a relatively 
large time-step (one tenth of the epicycle period) and plot the posi- 
tion of the particle at every time-step in Figure[T] As expected, SEI 
and SEKI follow the analytic solution exactly. The Quinn et al. in- 
tegrator yields an approximate ellipse but exhibits a phase error of 
6° per epicycle period. All three versions of the leapfrog integrator 
diverge badly from the analytic solution in less than one-quarter of 
an epicycle period. 

In Figure [2a] we plot the relative energy error as a function of 
the time-step. As expected, SEI and SEKI are exact to machine pre- 
cision for all time-steps. The Quinn et al. integrator and the mod- 
ified and symmetrized leapfrog integrators converge quadratically 
until machine precision is reached. Leapfrog converges only lin- 
early. In Figure |2b] we plot the computation time as a function of 
the relative error. This plot shows the fastest integrator for a desired 
precision. Of course, in this test case the choice is trivial, as SEI and 
SEKI give the exact solution for any time-step. 



4.3 Strongly perturbed epicyclic motion 

We also performed tests in which we place the test particle on an 
orbit with an impact parameter of IrHm, so the perturbing forces are 
much stronger than in the previous case. The trajectory of the test 
particle is then a horse-shoe orbit. 

In Figure|4]we plot the same diagnostics as in Figure|3]for this 
case. One can see that SEI and the Quinn et al. integrator perform 
equally well as measured by the energy error. However, the phase 
error is more than two orders of magnitude smaller using the SEI 
or the SEKI integrator rather than for any of the other integrators. 
Note that the SEKI integrator has significantly larger energy error 
than SEI. This is because for this test case all of the Hamiltonian 
operators in Eqs. and are roughly of the same magnitude 
and thus, the commutators of the operators that give rise to large 
integration errors. Because SEKI is based on a split into five oper- 
ators while SEI is based on three, there are more commutators and 
the total error is larger in SEKI. 

4.4 Perturbed Keplerian motion 

We finally test the strong gravity regime, where the forces that gov- 
ern the epicyclic motion can be viewed as a perturbation to a Keple- 
rian orbit. The initial positions and velocities are r = (0.125, 0, 0)^ 
and r = (0,-0.354,0)^, which correspond to an initially circular, 
bound orbit of m\ and mi- The semi-major axis of this orbit cor- 
responds to 0.18rHiii. The integration time is ten epicycle periods, 
corresponding to 226 Kepler orbital periods. 
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10*' 10"^ lO"'* 10"^ 10"^ 10^ 10° 



timestep [2it£2"^] 

(a) Relative energy error and phase error as a function of time-step. Lower 
is better 

Figure 3. Perturbed epicyclic motion (Sect. [42) . The trajectory passes the pe 
epicycle periods. 



phase error 

10"'° 10"^ 10° lO"* 10"^ 10° 




10 '° 10 ° 10 ° 10"" 10"^ 10° 

relative energy error 



(b) Computation time as a function of relative energy error and phase error. 
Lower is better. 

bing point mass at an impact parameter of 8 ruiu. The integration time is 100 



© 201 1 RAS, MNRAS 000,[T](9] 



Symplectic integrators in the shearing sheet 7 



phase error 

10'^° 10"^ 10*^ 10" 10"^ 10° 




lO '' 10"^ lO '* 10"^ 10 ^ 10 ' 10° 10 '° 10 ** 10 *^ 10"" 10"^ 10° 

timestep [2itn '] relative energy error 



(a) Relative energy error and phase error as a function of time-step. Lower (b) Computation time as a function of relative energy error and phase error, 
is better Lower is better. 

Figure 4. Strongly perturbed epicyclic motion (Sect. [4.3} . The trajectory has an impact parameter of 1 rniu corresponding to a horse-shoe orbit. The integration 
time is 100 epicycle periods. 
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Figure 6. Relative energy error using different integrators and a time-step 
At = 10"' IttQT^. The algorithm described in Appendix |a] was used to 
minimize round-off errors. 

In Figure|5]we plot again the same diagnostics as in Figure|2] 
The SEKI integrator is the most robust integrator in the sample, 
with the energy error at a given time-step more than two orders of 
magnitude smaller than its competitors. This performance comes 
with a drawback, as each time-step is computationally more expen- 
sive; nevertheless, at fixed energy error it is still more than one or- 
der of magnitude faster than its competitors. Its relative advantage 
improves further for orbits with smaller semi-major axis (relative 
to rHiu). 

Also note that for integrations in which the evaluation of 
forces is expensive (for example in a tree code), the SEKI integrator 
has a further advantage as it uses fewer time-steps. 

Figure |6] shows the time evolution of the relative energy error 
in the integration with At = 10"' 27rf2"' . All integrators (except the 
various flavors of leapfrog which are far off the scale and therefore 
not plotted) show no sign of a linear drift in the energy error: the 
maximum error is independent of time. This good behavior is due 
to the symplectic nature of these integrators as well as the use of 
the procedure described in Appendix |A] to control round-off error. 
The SEKI integrator is more than two orders of magnitude better 
than any other integrator in this example. 



5 CONCLUSIONS 

We have presented two new integrators for studying the small-scale 
dynamics of disks using Hill's approximation. 

The first is a simple symplectic and time-reversible integra- 
tor for Hill's equations of motion that we call symplectic epicy- 
cle integrator (SEI). In the absence of small-scale forces such as 
the self-gravitational forces of the disk particles, SEI solves Hill's 
equations of motion (epicyclic motion) exactly; otherwise the error 
over a fixed time integral scales as O(eAf-) when the small-scale 
forces are 0(e). Numerical tests using a variety of measures have 
shown that SEI always converges much faster than various flavors 
of leapfrog, and often much faster than the |Quinn et al.| ( [2(310) inte- 
grator. For small 6 (e ~ Sirynn/r)^ ~ 0.01) the phase error can be up 
to seven orders of magnitude smaller than the phase error from the 
Quinn et al. integrator at the same time-step (Figure|3^). The com- 
putational cost per time-step of SEI and the Quinn et al. integrator 



are equivalent. Although SEI is simple to code, a C implementation 
can be downloaded at |http : //sns ■ ias . edu/~rein/| 

The second integrator is also symplectic and time-reversible 
and is called the symplectic epicycle-Kepler integrator (SEKI). 
This integrator is useful in following bound two-body orbits in the 
sheared sheet and in this case can yield errors that are several orders 
of magnitude smaller than SEI. 

These integrators can be generalized in several ways. First, 
higher order integrators can be constructed by concatenating SEI 
and SEKI steps of varying lengths (e.g., |Yoshida||1993] >. Second, 
although the discussion in this paper has, for simplicity, focused on 
the three-body problem, SEI can be applied to the Af-body problem 
using shear periodic boundary conditions (e.g., |Richardson| 1994| l. 
The force f in Eq. l[TJ can also describe gas drag on particles, al- 
though in this case the dynamics is not described by a Hamiltonian 
so the advantage of a symplectic, time-reversible integrator is less 
clear. 

The SEKI integrator can also be generalized beyond Hill's ap- 
proximation. For example, consider a test particle orbiting in the 
gravitational field of a binary star with masses Mi and M2. The 
Hamiltonian can be written as a sum of two Keplerian Hamiltoni- 
ans minus the kinetic energy, 

H (r, p) = //Kep.Ml (r, P) + ^&p.M2 (r, P) - ^Kin (P) • (19) 

The first two terms can be solved exactly, and the last term is simply 
a drift. Thus, the analog of Eq. jl7[ l provides a second-order accu- 
rate, symplectic and time-reversible integrator that is exact when 
the test-particle motion is dominated by the gravitational field from 
either body. 
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We tested both implementations of the rotation operator, 
Eq. ijTTJ and Eq. (|ATJ. As expected, the implementation using shear 
operators shows no sign of a linear drift in energy (Figure |6). In 
contrast, for the test case presented in Figure|6] the straightforward 
implementation of Eq. produces a linear drift of about 10"' 
after only 100 epicycle periods. The additional computational cost 
of implementing the rotations using shear operators is negligible, 
and in long integrations this refinement can dramatically improve 
the accuracy. 



APPENDIX A: ROUND-OFF ERROR DUE TO 
ROTATIONS 

Floating-point operations on a computer are subject to rounding 
errors. The IEEE standard for floating-point arithmetic (IEEE 754) 
specifies the rounding algorithm and ensures that the rounding error 
for additions and multiplications is quasi-random. Thus, the round- 
ing error should grow with time as 0(e ■ (?/Af)''^), where e is the 
machine precision. 

Both the SEI and SEKI integrators use rotations, which are 
written in Eqs. and \13\ as a rotation matrix with angle (p = 
Q.At. However, operating with this matrix on a position vector (x, y) 
will not in general preserve the norm r = (x^+j^)''^ because numer- 
ically sin" ifi + cos' ifi t I. Since is the same at every time-step, the 
error grows linearly with the number of operations, 0{e ■ t), much 



worse than the 0{e ■ f''-) behavior described above ( Petit 1998 and 
references therein). This results in an undesirable linear drift in en- 
ergy. 

For other integrators that solve the Kepler problem using rota- 
tions, such as the Wisdom-Holman integrator, this problem is usu- 
ally not important. This is because the rotation angle in the Kepler 
problem depends implicitly on radius (Kepler's law) and thus the 
rotation is actually a so-called twist map. There exists a KAM-like 
theorem ( [Blank etal.|1997[ l for twist maps that restricts the solution 
to an invariant torus in phase space. 

|Petit| ( |1998^ describes one way to solve this problem by de- 
composing the rotation operator in Eq. ([TTJ into three shear opera- 
tors 



cos cp sin q 
- sin d> cos < 



anU 1 ' 1 ' -tani<4 1 



(Al) 



Why does this help? Suppose we replace the rotation matrix 
in Eq. l[TTJ by an arbitrary matrix R. It is then straightfor- 
ward to show that the transformation from (x" ,y" , p", p''.y to 
(x"+',y'+',pf' ,p';+')^ defined by Eqs. ((9| to iH^ is symplectic if 



and only if detR = 1. Since 1 and are represented exactly in 
floating-point arithmetic, each of the matrices on the right side of 
Eq. |ATJ has a determinant of exactly 1 . Thus the transformation is 
symplectic whether or not sin (/> and tan ^(f) are related by the appro- 
priate trigonometric identity, and hence is insensitive to round-off 
errors in evaluating these functions. 
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